clear all;

% First, compute the optimal reserve price for this distribution
m = .5;
v = .5^2;
mu = log((m^2)/sqrt(v+m^2));
sigma = sqrt(log(v/(m^2)+1));
% mu and sigma are thus the parameters of the lognormal distribution with mean m and variance v

fun = @(x)(-max(0.00001,x)*(1-normcdf(log(max(0.00001,x)),mu,sigma)));
x0 = 1;
rstar = fminsearch(fun,x0)
%rstar is the optimal reserve price for this distribution

%Entries in Table 2 are generated using Monte Carlo simulations

%First row

n = 2

nobs = 10000000; 
vals = lognrnd(mu,sigma,nobs,n);
vals = sort(vals,2);

alpha = .7;
refctr = [1];
for i = 1:n-1
    refctr = [refctr;alpha^i];
end;
    
aux = repmat(sort(refctr'),nobs,1);

refctr2 = refctr.*(1:n)';
refctr2 = (1-alpha)*refctr2(1:n-1);
aux2 = repmat(fliplr(refctr2'),nobs,1);

r=0;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
e0 = mean(sum(revsdirect,2)+aux3)

r=.1;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
e10 = mean(sum(revsdirect,2)+aux3)

r=(.1+rstar)/2;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
eMP = mean(sum(revsdirect,2)+aux3)

r=rstar;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
eOPT = mean(sum(revsdirect,2)+aux3)

%Second row

n = 6

nobs = 10000000; 
vals = lognrnd(mu,sigma,nobs,n);
vals = sort(vals,2);

alpha = .7;
refctr = [1];
for i = 1:n-1
    refctr = [refctr;alpha^i];
end;
    
aux = repmat(sort(refctr'),nobs,1);

refctr2 = refctr.*(1:n)';
refctr2 = (1-alpha)*refctr2(1:n-1);
aux2 = repmat(fliplr(refctr2'),nobs,1);

r=0;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
e0 = mean(sum(revsdirect,2)+aux3)

r=.1;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
e10 = mean(sum(revsdirect,2)+aux3)

r=(.1+rstar)/2;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
eMP = mean(sum(revsdirect,2)+aux3)

r=rstar;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
eOPT = mean(sum(revsdirect,2)+aux3)

%Third row

n = 10

nobs = 10000000; 
vals = lognrnd(mu,sigma,nobs,n);
vals = sort(vals,2);

alpha = .7;
refctr = [1];
for i = 1:n-1
    refctr = [refctr;alpha^i];
end;
    
aux = repmat(sort(refctr'),nobs,1);

refctr2 = refctr.*(1:n)';
refctr2 = (1-alpha)*refctr2(1:n-1);
aux2 = repmat(fliplr(refctr2'),nobs,1);

r=0;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
e0 = mean(sum(revsdirect,2)+aux3)

r=.1;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
e10 = mean(sum(revsdirect,2)+aux3)

r=(.1+rstar)/2;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
eMP = mean(sum(revsdirect,2)+aux3)

r=rstar;
gr = (vals>r);
cnt = sum(gr,2);
aux3 = r*(alpha.^(cnt-1)).*(cnt);
revsdirect = aux2.*vals(:,1:n-1).*gr(:,1:n-1);
eOPT = mean(sum(revsdirect,2)+aux3)
